# Clear memory
rm(list = ls())

# Load packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  maps, tidyverse, tidylog, data.table, janitor,
  broom, haven, sf, statar, fixest, patchwork, showtext, 
  rmapshaper, maps
)
options("tidylog.display" = NULL)

font_add_google("Lato")
showtext_auto()

########################################
## COUNTERFACTUAL MAPS
########################################

fips <- fread("data/simulation_fips_crosswalk.csv")

non_attainment <- fread("data/productivity_panel.csv") %>%
  filter(year >= 1986 & year <= 1997) %>%
  group_by(fips, industry) %>%
  mutate(total_nonattainment = sum(nonattainment, na.rm = T)) %>%
  ungroup() %>%
  mutate(nonattainment_1991 = nonattainment > 0 & year == 1991) %>%
  mutate(nonattainment_1997 = nonattainment > 0 & year == 1997) %>%
  group_by(fips) %>%
  mutate(nonattainment_1991 = max(nonattainment_1991)) %>%
  mutate(ever_nonattainment = max(nonattainment_1997) > 0) %>%
  ungroup() %>%
  filter(year == 1997) %>%
  mutate(fips = str_pad(fips, 5, "left", "0")) %>%
  distinct(fips, nonattainment_1997)

########################################
## LOAD PREREQS

# Load the county data from the maps package
state_fips <- state.fips %>%
  distinct(fips, .keep_all = TRUE)

dup_names <- str_sub(state_fips$polyname, 1, str_locate(state_fips$polyname, ":")[, 1] - 1)

state_fips$polyname[!is.na(dup_names)] <- dup_names[!is.na(dup_names)]

state_sdf <- st_as_sf(maps::map("state", plot = FALSE, fill = TRUE)) %>%
  rename(polyname = ID) %>%
  full_join(state_fips) %>%
  st_transform(5070)
state_sdf_simpler = rmapshaper::ms_simplify(state_sdf, keep_shapes = TRUE, keep = .05)

county_fips <- county.fips
counties_sdf <- st_as_sf(maps::map("county", plot = FALSE, fill = TRUE)) %>%
  rename(polyname = ID) %>%
  full_join(county_fips) %>%
  mutate(fips = ifelse(fips == 46113, 46102, fips)) %>%
  mutate(fips = ifelse(polyname == "florida,okaloosa", 12091, fips)) %>%
  mutate(fips = ifelse(polyname == "virginia,accomack", 51001, fips)) %>%
  mutate(fips = ifelse(polyname == "texas,galveston", 48167, fips)) %>%
  mutate(fips = ifelse(polyname == "louisiana,st martin", 22099, fips)) %>%
  mutate(fips = ifelse(polyname == "north carolina,currituck", 37053, fips)) %>%
  mutate(fips = ifelse(polyname == "washington,pierce", 53053, fips)) %>%
  mutate(fips = ifelse(polyname == "washington,san juan", 53055, fips)) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  st_transform(5070)

counties_sdf_simpler = rmapshaper::ms_simplify(counties_sdf, keep_shapes = TRUE, keep = .05)


########################################
## BASE RESULTS

simulation_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities = (exp(amenities - 1) - 1) * 100,
    output_change = nominal_wage * base_labor / cf_labor * 100,
    real_wage = (exp(real_wage - 1) - 1) * 100,
    emissions = (emissions_pm25 - 1) * 100,
    emissions = ifelse(industry == "other", -9999, emissions)
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(
    welfare = weighted.mean(welfare, base_labor),
    amenities = weighted.mean(amenities, base_labor),
    real_wage = weighted.mean(real_wage, base_labor),
    labor = sum(labor) / sum(base_labor) * 100,
    emissions = max(emissions)
  ) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070) |>
  st_simplify(dTolerance = 1e3)

############################
## total welfare
############################

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = welfare, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Welfare (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(-5, -2.5, 0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/welfare_map.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-4-welfare_map.eps"),
  width = 20,
  height = 10,
  dpi = 175
)

############################
## amenities welfare
############################

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = amenities, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Amenities\nWelfare (%)",
    colors = c("#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(0, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/welfare_amenities_map.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-5-bottom-left-welfare_amenities_map.eps"),
  width = 20,
  height = 10,
  dpi = 175
)
############################
## consumption welfare
############################

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = real_wage, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Consumption\nWelfare (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(-5, -2.5, 0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/welfare_consumption_map.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)

ggsave(
  paste0("output/figure-5-bottom-right-welfare_consumption_map.eps"),
  width = 20,
  height = 10,
  dpi = 175
)
############################
## population change
############################

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = labor, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Population (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(-5, -2.5, 0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/population_map.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-6-top-left-population_map.eps"),
  width = 20,
  height = 10,
  dpi = 175
)
############################
## labor reallocation welfare
############################

realloc_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities - amenities - 1,
    real_wage = real_wage - 1
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(welfare_realloc = weighted.mean(welfare, base_labor)) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)

simulation_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-off_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities - amenities - 1,
    real_wage = real_wage - 1
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(welfare_no_realloc = weighted.mean(welfare, base_labor)) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(realloc_results) %>%
  mutate(realloc_value = welfare_realloc - welfare_no_realloc)

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = realloc_value, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Labor\nReallocation\nWelfare (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(-5, -2.5, 0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )


ggsave(
  paste0("output/welfare_map_realloc_value.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-6-top-right-welfare_map_realloc_value.eps"),
  width = 20,
  height = 10,
  dpi = 175
)
realloc_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities - amenities - 1,
    real_wage = (exp(real_wage - 1) - 1) * 100
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(
    welfare_realloc = weighted.mean(welfare, base_labor),
    real_wage = weighted.mean(real_wage, base_labor)
  ) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)

simulation_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-off_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities - amenities - 1,
    real_wage = real_wage - 1
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(welfare_no_realloc = weighted.mean(welfare, base_labor)) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(realloc_results) %>%
  mutate(realloc_value = welfare_realloc - welfare_no_realloc)

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = realloc_value, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Trade\nReallocation\nWelfare (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-.04, .04),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-.04, -.02, -.01, 0, .01, .02, .04)),
    oob = scales::squish,
    breaks = scales::pretty_breaks(n = 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/welfare_map_trade_value.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-7-left-welfare_map_trade_value.eps"),
  width = 20,
  height = 10,
  dpi = 175
)
realloc_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities - amenities - 1,
    real_wage = real_wage - 1
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(welfare_realloc = weighted.mean(welfare, base_labor)) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)

simulation_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-off_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities - amenities - 1,
    real_wage = real_wage - 1
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(welfare_no_realloc = weighted.mean(welfare, base_labor)) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(realloc_results) %>%
  mutate(realloc_value = welfare_realloc - welfare_no_realloc)

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = realloc_value, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Pollution\nReallocation\nWelfare (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(-5, -2.5, 0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/welfare_map_physical_value.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-7-right-welfare_map_physical_value.eps"),
  width = 20,
  height = 10,
  dpi = 175
)


cong_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-on_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  group_by(fips) %>%
  summarise(
    base_labor = sum(base_labor),
    cf_labor = sum(cf_labor)
  ) %>%
  mutate(welfare = (exp((base_labor / cf_labor)^(-.3) - 1) - 1) * 100) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)

simulation_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities - amenities - 1,
    real_wage = real_wage - 1
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(welfare_no_cong = weighted.mean(welfare, base_labor)) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(cong_results) %>%
  mutate(cong_value = welfare - welfare_no_cong)

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070)


plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = welfare, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Welfare\nChange From\nCongestion and\nAgglomeration\n (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(-5, -2.5, 0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.98, 0.30)
  )

ggsave(
  paste0("output/welfare_map_congestion_value.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-e2-welfare_map_congestion_value.eps"),
  width = 20,
  height = 10,
  dpi = 175
)

############################
## manufacturing welfare
############################

simulation_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  filter(industry == "manufacturing") %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities = (exp(amenities - 1) - 1) * 100,
    output_change = nominal_wage * base_labor / cf_labor * 100,
    real_wage = (exp(real_wage - 1) - 1) * 100,
    emissions = (emissions_pm25 - 1) * 100,
    emissions = ifelse(industry == "other", -9999, emissions)
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(
    welfare = weighted.mean(welfare, base_labor),
    amenities = weighted.mean(amenities, base_labor),
    real_wage = weighted.mean(real_wage, base_labor),
    labor = sum(labor) / sum(base_labor) * 100,
    emissions = max(emissions)
  ) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)


plot_sdf <- counties_sdf_simpler %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070)


ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = welfare, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Manufacturing\nWelfare (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(-5, -2.5, 0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/welfare_map_manu.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-5-top-left-welfare_map_manu.eps"),
  width = 20,
  height = 10,
  dpi = 175
)
############################
## manufacturing population
############################

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = labor, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Manufacturing\nPopulation (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(-5, -2.5, 0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.96, 0.30)
  )

ggsave(
  paste0("output/manufacturing_pop_map.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)

ggsave(
  paste0("output/figure-6-bottom-left-manufacturing_pop_map.eps"),
  width = 20,
  height = 10,
  dpi = 175
)

############################
## nonmanufacturing welfare
############################

simulation_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  filter(industry == "other") %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities = (exp(amenities - 1) - 1) * 100,
    output_change = nominal_wage * base_labor / cf_labor * 100,
    real_wage = (exp(real_wage - 1) - 1) * 100,
    emissions = (emissions_pm25 - 1) * 100,
    emissions = ifelse(industry == "other", -9999, emissions)
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(
    welfare = weighted.mean(welfare, base_labor),
    amenities = weighted.mean(amenities, base_labor),
    real_wage = weighted.mean(real_wage, base_labor),
    labor = sum(labor) / sum(base_labor) * 100,
    emissions = max(emissions)
  ) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)

plot_sdf <- counties_sdf_simpler %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = welfare, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Nonmanufacturing\nWelfare (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(-5, -2.5, 0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.99, 0.30)
  )

ggsave(
  paste0("output/welfare_map_nonmanu.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)

ggsave(
  paste0("output/figure-5-top-right-welfare_map_nonmanu.eps"),
  width = 20,
  height = 10,
  dpi = 175
)


############################
## nonmanufacturing population
############################

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = labor, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Nonmanufacturing\nPopulation (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = c(-5, -2.5, 0, 2.5, 5)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.98, 0.30)
  )

ggsave(
  paste0("output/nonmanufacturing_pop_map.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)

ggsave(
  paste0("output/figure-6-bottom-right-nonmanufacturing_pop_map.eps"),
  width = 20,
  height = 10,
  dpi = 175
)


############################
## optimal vs actual
############################

simulation_results <- read.csv("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(1 / welfare - 1) - 1) * 100,
    amenities = (exp(amenities - 1) - 1) * 100,
    output_change = nominal_wage * base_labor / cf_labor * 100,
    real_wage = (exp(real_wage - 1) - 1) * 100,
    emissions = (emissions_pm25 - 1) * 100,
    emissions = ifelse(industry == "other", -9999, emissions)
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(
    welfare = weighted.mean(welfare, base_labor),
    amenities = weighted.mean(amenities, base_labor),
    real_wage = weighted.mean(real_wage, base_labor),
    labor = sum(labor) / sum(base_labor) * 100,
    labor_level = sum(base_labor),
    emissions = max(emissions)
  ) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)


plot_sdf <- counties_sdf_simpler %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070) |>
  mutate(worse_off = welfare < 0)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = welfare, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "First-Best\nVersus\nActual (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-5, 5),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-5, -2, -1, 0, 1, 2, 5)),
    oob = scales::squish,
    breaks = scales::pretty_breaks()
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.99, 0.30)
  )

ggsave(
  paste0("output/welfare_map_opt.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-e4-welfare_map_opt.eps"),
  width = 20,
  height = 10,
  dpi = 175
)
############################
## productivity effects
############################

cong_results <- read.csv("data/counterfactual/welfare_opt-base_prod--0.03_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities = (exp(amenities - 1) - 1) * 100,
    output_change = nominal_wage * base_labor / cf_labor * 100,
    real_wage = (exp(real_wage - 1) - 1) * 100,
    emissions = (emissions_pm25 - 1) * 100,
    emissions = ifelse(industry == "other", -9999, emissions)
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(prod_welfare = weighted.mean(welfare, base_labor)) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)

simulation_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities = (exp(amenities - 1) - 1) * 100,
    output_change = nominal_wage * base_labor / cf_labor * 100,
    real_wage = (exp(real_wage - 1) - 1) * 100,
    emissions = (emissions_pm25 - 1) * 100,
    emissions = ifelse(industry == "other", -9999, emissions)
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(welfare = weighted.mean(welfare, base_labor)) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) |>
  left_join(cong_results) %>%
  mutate(prod_value = prod_welfare - welfare)

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = prod_value, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Welfare Change\nfrom -3%\nProductivity\nShock (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-.25, .25),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-.25, -.2, -.05, 0, .05, .1, .25)),
    oob = scales::squish,
    breaks = c(-.2, -.1, 0, .1, .2)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.98, 0.30)
  )

ggsave(
  paste0("output/welfare_map_prod_-3.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-e3-left-welfare_map_prod_-3.eps"),
  width = 20,
  height = 10,
  dpi = 175
)
cong_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.03_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities = (exp(amenities - 1) - 1) * 100,
    output_change = nominal_wage * base_labor / cf_labor * 100,
    real_wage = (exp(real_wage - 1) - 1) * 100,
    emissions = (emissions_pm25 - 1) * 100,
    emissions = ifelse(industry == "other", -9999, emissions)
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(prod_welfare = weighted.mean(welfare, base_labor)) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  left_join(non_attainment)

simulation_results <- read.csv("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  mutate(
    welfare = (exp(welfare - 1) - 1) * 100,
    amenities = (exp(amenities - 1) - 1) * 100,
    output_change = nominal_wage * base_labor / cf_labor * 100,
    real_wage = (exp(real_wage - 1) - 1) * 100,
    emissions = (emissions_pm25 - 1) * 100,
    emissions = ifelse(industry == "other", -9999, emissions)
  ) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  group_by(fips) %>%
  summarise(welfare = weighted.mean(welfare, base_labor)) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) |>
  left_join(cong_results) %>%
  mutate(prod_value = prod_welfare - welfare)

plot_sdf <- counties_sdf_simpler %>%
  left_join(simulation_results) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  st_transform(5070)

ggplot() +
  geom_sf(
    data = plot_sdf,
    aes(fill = prod_value, color = nonattainment_1997, size = nonattainment_1997)
  ) +
  geom_sf(
    data = state_sdf_simpler,
    col = "#555555",
    size = 0.4,
    fill = "transparent"
  ) +
  scale_color_manual(values = c("#AAAAAA", "#000000"), guide = "none") +
  scale_size_manual(values = c(.1, .5), guide = "none") +
  scale_fill_gradientn(
    name = "Welfare Change\nfrom +3%\nProductivity\nShock (%)",
    colors = c("#b2182b", "#f4a582", "#fddbc7", "#ffffff", "#92c5de", "#4393c3", "#2166ac"),
    na.value = "#AAAAAA",
    aesthetics = "fill",
    limits = c(-.25, .25),
    guide = guide_colourbar(
      ticks.linewidth = 2,
      ticks.colour = "#333333",
      frame.colour = "#333333",
      frame.linewidth = 2
    ),
    values = scales::rescale(c(-.25, -.2, -.05, 0, .05, .1, .25)),
    oob = scales::squish,
    breaks = c(-.2, -.1, 0, .1, .2)
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    text = element_text(family = "Lato"),
    legend.title = element_text(size = 30),
    legend.text = element_text(size = 30),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1.5, "cm"),
    legend.position = c(0.98, 0.30)
  )

ggsave(
  paste0("output/welfare_map_prod_3.pdf"),
  width = 20,
  height = 10,
  dpi = 175
)
ggsave(
  paste0("output/figure-e3-right-welfare_map_prod_3.eps"),
  width = 20,
  height = 10,
  dpi = 175
)

############################
## ramp up nonattainment
############################


files <- list.files("data/counterfactual/", pattern = "welfare_total_heterog-off_perc*", full.names = T)[seq(1, 101, 5)]
welfare <- map_dfr(
  files,
  function(file) {
    temp <- read_csv(file) %>%
      mutate(percent = str_match(file, "percent-\\s*(.*?)\\s*.csv")[2])
  }
) %>%
  mutate(percent = as.numeric(percent)) %>%
  drop_na(percent) %>%
  mutate(
    max = total == max(total),
    across(total:consumption_total, ~ (exp(.x - 1) - 1) * 100)
  )

# set up rescaling for secondary axis
xlim.prim <- c(min(welfare$percent), max(welfare$percent))
xlim.sec <- c(min(welfare$total_in_na), max(welfare$total_in_na))
b <- diff(xlim.prim) / diff(xlim.sec)
a <- xlim.prim[1] - b * xlim.sec[1]

ylim.prim <- c(min(welfare$total), max(welfare$amenity_total))
ylim.sec <- c(min(welfare$consumption_total), max(welfare$consumption_total))
by <- diff(ylim.prim) / diff(ylim.sec)
ay <- ylim.prim[1] - by * ylim.sec[1]


welfare %>%
  ggplot() +
  geom_line(aes(x = percent * 100, y = total), size = 1.5, color = "black") +
  geom_point(aes(x = percent * 100, y = total), size = 6) +
  geom_line(aes(x = percent * 100, y = amenity_total), size = 1.5, color = "#2166ac", linetype = "dashed") +
  geom_point(aes(x = percent * 100, y = amenity_total), color = "#2166ac", size = 6) +
  geom_line(aes(x = percent * 100, y = ay + consumption_total * by), size = 1.5, color = "#b2182b", linetype = "dotdash") +
  geom_point(aes(x = percent * 100, y = ay + consumption_total * by), size = 6, color = "#b2182b") +
  scale_color_manual(values = c("black", "black")) +
  scale_shape(solid = T) +
  annotate(
    geom = "text", x = 92, y = .84, label = "Consumption", hjust = 0,
    color = "#b2182b", size = 10
  ) +
  annotate(
    geom = "text", x = 52, y = .85, label = "Amenities", hjust = 0,
    color = "#2166ac", size = 10
  ) +
  annotate(
    geom = "text", x = 30, y = .72, label = "Total", hjust = 0,
    color = "black", size = 10
  ) +
  # annotate(
  #     geom = "text", x = 33, y = .50, label = "Second-Best Threshold", hjust = 0,
  #     color = "#c51b7d", size = 10
  # ) +
  # annotate(
  #     geom = "segment", x = 15, xend = 15, y = .5225-.03, yend = welfare$total[welfare$max==T] + .007,
  #     color = "#c51b7d", size = 1
  # ) +
  annotate(
    geom = "text", x = 60, y = .775, label = "More Stringent Policy", hjust = 0,
    color = "black", size = 10
  ) +
  geom_segment(
    aes(x = 60, xend = 25, y = .76, yend = 0.76),
    size = 1.5, arrow = arrow(length = unit(0.15, "inches"))
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    title = element_text(size = 24),
    axis.text.x = element_text(size = 30), axis.text.y = element_text(size = 30),
    axis.title.x = element_text(size = 30), axis.title.y = element_text(size = 30),
    panel.grid.minor.x = element_blank(), panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank(),
    axis.line = element_line(colour = "black"), axis.ticks = element_line(),
    axis.title.y.right = element_text(vjust = 2), axis.title.y.left = element_text(vjust = 2),
    axis.title.x.top = element_text(vjust = 1), plot.margin = margin(15, 15, 15, 15)
  ) +
  labs(
    y = "Total Welfare and Amenity Welfare (%)",
    x = "Counterfactual Nonattainment Threshold (% of 1997 Thresholds)"
  ) +
  scale_x_reverse(
    breaks = scales::pretty_breaks(),
    sec.axis = sec_axis(~ rev(. / 100 - a) / b,
      breaks = scales::pretty_breaks(),
      name = "Number of Counties in Nonattainment"
    )
  ) +
  scale_y_continuous(
    breaks = scales::pretty_breaks(n = 7),
    sec.axis = sec_axis(~ (. - ay) / by,
      name = "Consumption Welfare (%)",
      breaks = scales::pretty_breaks(n = 7)
    )
  )


ggsave("output/nonattainment_results.pdf", width = 15, height = 15)
ggsave("output/nonattainment_results.eps", width = 15, height = 15)

welfare %>%
  ggplot() +
  geom_line(aes(x = percent * 100, y = total), size = 1.5, color = "black") +
  geom_point(aes(x = percent * 100, y = total), size = 6) +
  geom_line(aes(x = percent * 100, y = amenity_total), size = 1.5, color = "black", linetype = "dashed") +
  geom_point(aes(x = percent * 100, y = amenity_total), color = "black", size = 6) +
  geom_line(aes(x = percent * 100, y = ay + consumption_total * by), size = 1.5, color = "black", linetype = "dotdash") +
  geom_point(aes(x = percent * 100, y = ay + consumption_total * by), size = 6, color = "black") +
  scale_color_manual(values = c("black", "black")) +
  scale_shape(solid = T) +
  annotate(
    geom = "text", x = 92, y = .84, label = "Consumption", hjust = 0,
    color = "black", size = 10
  ) +
  annotate(
    geom = "text", x = 52, y = .85, label = "Amenities", hjust = 0,
    color = "black", size = 10
  ) +
  annotate(
    geom = "text", x = 30, y = .72, label = "Total", hjust = 0,
    color = "black", size = 10
  ) +
  # annotate(
  #     geom = "text", x = 33, y = .50, label = "Second-Best Threshold", hjust = 0,
  #     color = "#c51b7d", size = 10
  # ) +
  # annotate(
  #     geom = "segment", x = 15, xend = 15, y = .5225-.03, yend = welfare$total[welfare$max==T] + .007,
  #     color = "#c51b7d", size = 1
  # ) +
  annotate(
    geom = "text", x = 60, y = .775, label = "More Stringent Policy", hjust = 0,
    color = "black", size = 10
  ) +
  geom_segment(
    aes(x = 60, xend = 25, y = .76, yend = 0.76),
    size = 1.5, arrow = arrow(length = unit(0.15, "inches"))
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    title = element_text(size = 24),
    axis.text.x = element_text(size = 30), axis.text.y = element_text(size = 30),
    axis.title.x = element_text(size = 30), axis.title.y = element_text(size = 30),
    panel.grid.minor.x = element_blank(), panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank(),
    axis.line = element_line(colour = "black"), axis.ticks = element_line(),
    axis.title.y.right = element_text(vjust = 2), axis.title.y.left = element_text(vjust = 2),
    axis.title.x.top = element_text(vjust = 1), plot.margin = margin(15, 15, 15, 15)
  ) +
  labs(
    y = "Total Welfare and Amenity Welfare (%)",
    x = "Counterfactual Nonattainment Threshold (% of 1997 Thresholds)"
  ) +
  scale_x_reverse(
    breaks = scales::pretty_breaks(),
    sec.axis = sec_axis(~ rev(. / 100 - a) / b,
                        breaks = scales::pretty_breaks(),
                        name = "Number of Counties in Nonattainment"
    )
  ) +
  scale_y_continuous(
    breaks = scales::pretty_breaks(n = 7),
    sec.axis = sec_axis(~ (. - ay) / by,
                        name = "Consumption Welfare (%)",
                        breaks = scales::pretty_breaks(n = 7)
    )
  )

ggsave("output/figure-3-nonattainment_results_bw.eps", width = 15, height = 15)


